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(57) Abstract: The instant invention teaches an improved frequency domain based method of generating attributes trom 2D or 3D 
seismic data. The central idea of the instant invention is that if a curve characterized by constant coefficients is filled to a short- 
window Fourier transform amplitude or phase spectrum, the coefficient estimates that arc obtained thereby can be used as seismic 
attributes for use in detecting changes in subsurface properties that may be associated with the generation, migration, accumulation, 

^ or presence of hydrocarbons. Since it is well known that many subsurface structural and stratigraphic features alter the frequency 
characteristics of seismic waveforms passing therethrough, the instant invention provides a new way to search through large numbers 

1^ of seismic traces for frequency changes that may be correlated with subsurface features of exploration interest. 
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Method of Frequency Domain 

Seismic Attribute Generation 
TECHNICAL FIELD 

The present invention relates genially to tiie field of seismic exploration and, in more 
particular, to methods of quantifying and visualizing subsurfece structural and stratigraphic 
features in two and three dimensions. TTiis invention also relates to the field of seismic 
attribute generation and the use of seismic attributes to detect tfie presence of hydrocarbons in 
the subsurface. Additionally, it relates to the correlation of seismic attributes with subsurface 
features that are conducive to the migration, accumulation, and presence of hydrocarbons. 
The mvention disclosed herein will be most fully appreciated by those in the seismic 
interpretation and seismic processing arts. 

BACKGROUND 

A seismic survey represents an attempt to map the subsurface of the earth by sending 
sound energy down into tiie ground and recording the "echoes" that return from the rock 
layers below. The source of the down-going sound energy might come, for example, from 
explosions or seismic vibrators on land, or air guns in marine environments. During a 
seismic survey, the energy source is moved to various positions across the surface of the earth 
above a geological structure of interest. Each time the source is activated, it generates a 
seismic signal that travels downward through the earth, is reflected, and, upon its return, is 
recorded at a great many locations on the surface. Multiple source / recording combinations 
are then combined to create a near continuous profile of the subsurface that can extend for 
many miles. In a two-dimensional (2D) seismic survey, the recording locations are generally 
laid out along a single straight line, whereas in a three dimensional (3D) survey the recording 
locations are distributed across the surface in a grid pattern. In simplest terms, a 2D seismic 
line can be thought of as giving a cross sectional picture (vertical slice) of the earth layers as 
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they exist directly beneath the recording locations. A 3D survey produces a data "cube" or 
volume that is, at least conceptually, a 3D picture of the subsurface that lies beneath the 
survey area. ]n reality, though, both methods interrogate some volume of the earth lying 
beneath the area covered by the survey. 

5 A seismic survey is composed of a very large number of individual seismic recordings or 

traces, hi a typical 2D survey, there will usually be several tens of thousands of traces, 
whereas in a 3D survey the number of individual traces may run into the multiple millions of 
traces. The term "unstacked" seismic traces is used by those skilled in the art to describe 
seismic traces as they are collected in field recordings. This t^m also is applied to seismic 

10 traces during the processing sequence up to the point where traces are "stacked" or averaged 
together. General background information pertaining 3D data acquisition and processing may 
be found in Chapter 6, pages 384-427, of Seismic Data Processing by Ozdogan Yilmaz, 
Society of Exploration Geophysicists, 1987, the disclosure of which is incorporated herein by 
reference. Chapter 1, pages 9 to 89, of Yilmaz contains general information relating to 

15 conventional 2D processing and that disclosure is also incorporated herein by reference. 

A modem seismic trace is a digital recording (analog recordings were used m the past) of 
the energy reflecting back from inhomogeneities or discontinuities in the subsurface, a partial 
reflection occurring each time there is a change in the elastic properties of the subsurfece 
materials. The digital samples are usually acquired at 0.002 second (2 millisecond or "ms") 

20 intervals, although 4 millisecond and 1 millisecond sampling intervals are also common. 
Thus, each digital sample in a seismic trace is associated with a travel time (in the case of 
reflected energy, a two-way travel time from the surface to the reflector and back to the 
surface again). Further, the surface location of each trace in a seismic survey is carefully 
recorded and remains associated with that trace during subsequent processing. This allows 
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the seismic information contained within the traces to be later correlated with specific surface 
and subsurfece locations, thereby providing a means for posting and contouring seismic data 
— and attributes extracted therefiom — on a map (i.e., "mapping"). 

The data in a 3D survey are amenable to viewing in a number of different ways. Firet, 
5 horizontal "constant time slices" may be extracted fi-om a stacked or unstacked seismic 
volume by collecting all of the digital samples that occur at the same travel time. This 
operation results in a 2D plane of seismic data. Sunilarly, a vertical plane of seismic data 
may be taken at an arbitrary azimuth through the volume by collecting and displaying the 
seismic traces that lie along a particular line. This operation, in efifect. extracts an individual 
10 2D seismic line fix)m within the 3D data volume. 

Seismic data that have been properly acquired and processed can provide a wealth of 
infomiation to the explorationist. one of the individuals within an oil company whose job it is 
to locate potential drilling sites. For example, a seismic profile gives the explorationist a 
broad view of the subsurfece structure of the rock layere and often reveals important features 
15 associated with the entrapment and storage of hydrocarbons such as faults, folds, anticlines, 
unconformities, and sub-suifece salt domes and reefe, among many others. During the 
computer processing of seismic data, estimates of subsurface rock velocities are routinely 
generated and near surfece inhomogeneities are detected and displayed. In some cases, 
seismic data can be used to directly estimate rock porosity, water saturation, and hydrocarbon 
20 content. Less obviously, seismic waveform attributes such as phase, peak amplitude, peak-to- 
trough ratio, and a host of others, can often be empirically correlated with known 
hydrocarbon occurrences and that correlation subsequently applied to seismic data collected 
over other exploration targets! 
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That being said, one of the most challenging tasks facing the seismic interpreter — 
one of the individuals within an oil company that is responsible for reviewing and analyzing 
the collected seismic data — is locating these stratigraphic and structural features of interest 
within individual seismic lines or, more problematically, within potentially enormous seismic 

5 volumes. By way of example only, it can be important for exploration purposes to locate 
regions in the subsurface in which the frequency content of seismic events reflected therefrom 
and transmitted therethrough are different from the surrounding rocks, as these oft-times 
subtle frequency changes may be indicative of the presence of fluids (including gas) in the 
subsurface rocks. Additionally, rock stratigraphic information may be revealed through the 

10 analysis of spatial variations in a seismic reflector^s character because these variations may 
often be empirically correlated with changes in reservoir lithology or fluid content. Since the 
precise physical mechanism which gives rise to these variations may not be well understood, 
it is common practice for interpreters to calculate a variety of attributes from the recorded 
seismic data and then plot or map them, looking for an attribute that has some predictive 

15 value. Given the enormous amount of data collected in a 3-D volume, further automated 
methods of enhancing the appearance of subsurface features related to the migration, 
accumulation, and presence of hydrocarbons are sorely needed. 

Before proceeding to a description of the present invention, however, it should be noted 
and remembered that the description of the invention which follows, together with the 

20 accompanying drawings, should not be construed as limiting the invention to the examples 
(or preferred embodiments) shown and described. This is so because those skilled in the art 
to which the invention pertains will be able to devise other forms of this invention within the 
ambit of the appended claims. 
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SUMMARY OF THE INVENTION 

The instant inventor has discovered an improved frequency domain based method of 
generating attributes from 2D or 3D seismic data. The central idea of the instant invention is 
that if a curve characterized by one or more constant coefficients is fitted to a Fourier 
transform amplitude or phase spectrum that has been calculated over a relatively short time 
window, the coefiBcient estimates that are obtained thereby are seismic attributes that can be 
used to detect changes in subsurface properties that may be associated with the generation, 
migration, accumulation, or presence of hydrocarbons. Since it is well known that many 
subsurfece structural and stratigraphic features alter the frequency characteristics of seismic 
waveforms passing therethrough, the instant invention provides a new way to search through 
large numbers of seismic traces for frequency changes that may be correlated with subsurfece 
features of exploration interest. 

As a preluninary step, the instant mefliod typically begins with the selection of an 
input 2D or 3D seismic data set that has been collected over a predetermined volume of the 
earth that contains structural or stratigraphic features of mterest to the explorationist. The 
selected seismic traces may either be stacked or unstacked, however stacked seismic traces 
are the preferred input data for use with this invention. 

As a next step, a zone of interest within seismic data is specified and an initial 
analysis window within the zone of interest is selected. The zone of interest is preferably 
defined in terms of time (or depth) and lateral extent within the selected seismic data set The 
time span of the analysis window will typically be somewhat smaller than that of the zone of 
interest and the instant inventor contemplates that a series of sliding analysis windows will be 
used to temporally "cover" the entire zone of interest. In the case of a 3D data set, each 
analysis window ultimately gives rise to a 2-D "plane" of coefficients — one coefficient 
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being supplied by each seismic trace in the volume — so that, by performing the steps that 
follow on a large number of windows, a 3-D volume may be constructed by combining the 
resulting 2-D planes. When the instant method is applied to 2-D seismic data, the ou^ut data 
set (assuming that multiple analysis windows are employed) is again a 2-D data set. 

5 Within each analysis window, a Fourier transform phase or amplitude spectrum is 

next calculated using the samples contained therein. The resulting spectrum is then fit by a 
fimction characterized by constant coefficients, thereby resulting in a collection of coefficient 
estimates that are associated with this spectrum and function. One of the coefficient 
estimates from the curve fit is then selected and displayed at the same relative location (in 

10 space and time) as the center of the analysis window. This coefficient estimate reflects 
(depending on which particular coefficient is chosen) some feature related to die overall 
shape of the spectrum and, hence, the frequency content of the seismic trace fi*om which it 
was calculated. For example, fi-equency spectra that are front (i.e., low fi-equency ) loaded 
can be differentiated from those that are rear (i.e., high frequency) loaded because those 

15 spectra have different overall shapes, which shapes will be reflected in the estimated 
coefficient values. Further, spectra that are ununodal in shape can be differentiated from 
those that are bimodal (e.g., those spectra that contain a pronounced frequency notch). 
Clearly, many other variations are possible. 

If a number of seismic traces / analysis windows are processed by the previous 

20 method, the resulting coefficients can be assembled into seismic traces, lines, slices / planes 
and volumes which are suitable for viewing by the explorationist. Depending on the 
particular fitting function and coefficient selected for viewing, the explorationist will in a 
position to quantitatively assess changing frequency content and changes spectral shape, and 
to correlate those changes with the subsurface lithology. Since it is well Icnown that 
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subsurface geology can induce changes in the fi?equency content of seismic data, this method 
provides a way for the explorationist to rapidly scan large volumes of seismic data for these 
sorts of changes. 

Further, rather than displaying just individual coefiBcient estimates, mathematical 
5 combinations of any number of the calculated coefBcients at the same time level and trace 
may be formed and displayed. Still further, mathematical operations may be performed using 
coefficient estimates from different time levels in the same trace or between the same time 
levels on adjacent traces. 

Still further, any statistic related to the quality of the curve fit (e.g., statistical 
10 correlation or "R^") can also be calculated and used as a seismic attribute. Even forther, the 
display of the resulting seismic attribute might be conditioned on how well each spectrum is 
fit by the selected equation with, for example, atti-ibutes correspondmg to "good" fits being 
displayed while those corresponding to "bad" fits being set to some null value. This would 
help the explorationist judge whether a particular frequency change was truly representative 
15 ofthe underlying data or was a spurious result based on noise. 

Each of ti»e coefficients that results from tiie above fitting process is a seismic 
attribute that may be displayed and empirically or automatically correlated with subsurface 
structure or sti-atigraphy. The output from the instant metiiod is a new seismic atfribute that 
quantifies tiie shape ofthe spectrum and yields an atti-ibute tiiat numerically characterizes tiiat 
20 shape, hi the past, these frequency spectra could only be examined one-at-a-time, if the goal 
was to ascertain the overall shape of each spectrum. However, the instant metiiod provides a 
way to look at large numbers of spectral shapes in their proper spatial relationship and to 
delimit the areal extent of particular shapes. This method can be used on stacked or 
unstacked seismic data, 2D or 3D. 
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While this invention is susceptible of embodiment in many different forms, there is 
shown in the drawings, and will hereinafter be described in detail, one or more specific 
embodiments of the instant invention. It should be understood, however, that the present 
disclosure is to be considered an exemplification of the principles of the invention and is not 
5 intended to limit the invention to any specific embodiments or algorithms so described. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 illustrates the general environment in which the invention disclosed herein 
10 would typically be used. 

Figure 2 is a schematic diagram that illustrates how the instant invention might be 
utilized within a conventional exploration seismic processing stream. 

Figure 3 contains examples of whole-trace and short-window Fourier transforms. 
Figure 4 is a schematic diagram that illustrates how the instant invention operates to 
15 produce a plane of seismic attributes. 

Figure 5 illustrates how planes of seismic attributes produced by the instant method 
may be combined to form a seismic volume. 

DETAILED DESCRIPTION 

20 

Environment of the Invention 
Figure 1 illustrates the general environment in which the instant invention would 
typically be used. Seismic data 110 (either 2-D or 3-D) are collected in the field over a 
subsurface target of potential economic importance and are then taken to a processing center, 
25 where a variety of preparatory processes 120 are applied to the seismic traces to make them 
ready for use by the methods disclosed hereinafter. The preparatory processes typically 
include the association of an X and Y surface coordinate with every processed trace (e.g., by a 
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geometry routine). The processed traces 130 would then be made available for use by the 
instant invention and might be stored, by way of example only, on hard disk, magnetic tape, 
magneto-optical disk, DVD disl^ or other mass storage means. 

The methods disclosed herein would best be implemented in the form of a compiled 
5 computer program 140 that has been loaded onto a general purpose programmable computer 
150 where it is accessible by a seismic interpreter or processor. For purposes of the instant 
disclosure, a general purpose computer 150 will be defined to include, in addition to 
mainframes and woricstations, computea-s Aat provide for parallel or massively parallel 
computations, wherein the computational load is distributed between two or more processors. 
10 As is also illustrated in Figure 1, a digitized zone of interest is preferably also specified and 
provided as input 160 to the software that implements the instant methods. The exact means 
by which a zone of interest is determined, digitized, stored, and later read during program 
execution is unimportant to the instant invention and those ddlled in tiie art will recognize 
that this might be done any number of ways. 
15 A computer program 140 embodying the instant invention might be conveyed into the 

computer fliat is to execute it by means o^ for example, a floppy disk, a magnetic disk, a 
magnetic tape, a magneto-optical disk, an optical disk, a CD-ROM, a DVD disk, or loaded 
over a network. After the instant methods have been applied, the resulting seismic attribute 
traces would then typically be displayed either at a high resolution color computer monitor 
20 170 or in hard-copy form as a printed seismic section or a map 180. The traces feat result 
might be a 2D seismic attribute line, a "plane" of seismic attributes, or a seismic attribute 
volume, although the precise form the ou^ut data takes is not a part of tiie instant invention. 
The seismic interpreter would then use the displayed images to assist him or her in identifying 
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subsurface features conducive to the generation, migration, accumulation, or presence of 
hydrocarbons. 

Preparatory Processing 
As a first step, and as is generally illustrated in Figure 2, a 2D or 3D seismic survey is 

5 conducted over a particular subsurface earth volume (step 210). The data that are collected 
consist of unstacked (i.e., unsummed) seismic traces which contain digital information 
representative of rock layers within the volume of the earth lying beneath the survey. 
Methods by which such data are obtained in a form suitable for use by seismic processors and 
interpreters are well known to those skilled in the art. Additionally, those skilled in the art 

10 will recognize that the processing steps illustrated in Figure 2 are only broadly representative 
of the sorts of steps that seismic data would normally go through before it is interpreted: the 
choice and order of the processing steps, and the particular algorithms involved, may vary 
markedly depending on the particular seismic processor, the signal source (dynamite, 

I 

vibrator, etc.), the survey location (land, sea, etc.) of the data, and the company that processes 
15 the data. 

As is broadly illustrated in Figure 2, seismic data typically pass through a number of 
processing steps before it is ready to be used by interpreters. A common early step is the 
specification of the geometry of the survey (step 220), at which time each seismic trace is 
associated with both the position of the physical receiver (or array) on or near the surface of 
20 the earth that recorded that particular trace, and the "shof ' (or generated seismic signal) that 
provided the source for the recorded seismic energy. The shot-receiver location information 
is later used to determine the assigned position of the "stacked" seismic traces and for many 
other purposes. 
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After the initial pre-stack processing is completed, it is customary to condition the 
seismic signal on the unstacked seismic traces before creating a stacked (or summed) data 
volume (step 230). In Figure 2, the "Signal Processing / Conditioning / Imaging" step 230 
suggests a typical processing sequence, although those skilled in the art will recognize that 
many alternative processes could be used in place of the ones listed in the figure. 

An ultimate goal of seismic processmg is usually to produce a stacked seismic volume 
240 or a stacked seismic line (2-D data). As is well known to those skilled in the art. a CMP 
gather is a collection of traces that all have the same midpoint between a shot and a receiver. 
A stack (or CMP stack) is horizontal sum of all the traces in a CMP gather, i.e., the numerical 
values from each trace at the same time point are summed together to produce a single output 
"stacked" trace that is an "average" of the input traces. Additionally, before stacking the 
traces it is customary to first correct them for offeet from the shot ("normal moveouf or 
NMO. step 220). however the NMO correction will be assumed herein to be incorporated as 
part of the CMP stack and will not be considered separately. 

The explorationist may do an initial interpretation 250 of the resulting stacked 
volume, wherein he or she locates and identifies the principal reflectors and faults wherever 
they occur in the data set Finally, as noted in Figure 2. the explorationist will then typically 
use the processed seismic data to locate subsurface sfructural or stratigraphic features 
conducive to the generation, accumulation, or migration of hydrocarbons (i.e., prospect 
generation 290). This effort may incorporate additional data from a variety of non-seismic 
sources including, for example, well logs, satellite surveys, gravity surveys, etc. 

Another important use for seismic data is as a source for seismic attributes (step 270). 
As is well known to those skilled in the art, seismic attributes are values that are calculated 
from the seismic data and that serve to highlight some specific property or feature of the data 
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that might not otherwise be observable within the original data. Although Figure 2 seems to 
indicate that seismic attribute generation 270 takes place relatively late in the processing 
sequence, that is not always the case and attributes might potentially be calculated at almost 
any point in Figure 2. Step 270 is one point in a typical seismic processing sequence that 

5 could include the methods of the instant disclosure. Since seismic attributes can reveal 
subsurface details that are at odds with the initial seismic data interpretation 250, the 
explorationist may use the result of step 270 to help refine or reinterpret (step 280) the 
original interpretation (250) before moving on to the prospect generation 290 stage. 
As is suggested in Figure 2, any digital sample within a seismic volume 240 is 

10 uniquely identified by an (X,Y,TIME) triplet: the X and Y coordinates representing some 
location on the surface of the earth, and the time coordinate measuring a distance down the 
seismic trace. For purposes of specificity only, it will be assumed that the X direction 
corresponds to the "in-line" direction, and the Y measurement corresponds to the "cross-line" 
direction, as the terms "in-line" and "cross-line" are generally understood to mean in the art. 

15 Although time is the preferred and most commonly encountered vertical axis unit, those 
skilled in the art understand that other units are certainly possible and might include, for 
example, depth or fi-equency. Additionally, it is well known to those skilled in the art that it 
is possible to convert seismic traces from one axis unit (e.g., time) to another (e.g., depth) 
using standard mathematical conversion techniques. 

20 For purposes of the instant invention, data that are suitable for analysis by the methods 

disclosed herein might consist of, for purposes of illustration only, an unstacked 2-D seismic 
line, an unstacked 2-D stacked seismic line extracted from a 3-D seismic survey , an 
unstacked 3-D portion of a 3-D seismic survey, or, preferably, a portion of a stacked 3-D 
survey. The invention disclosed herein is most effective when applied to a group of 
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unstacked seismic traces that have an unda-l>ing spatial relationship with respect to some 
subsurface geological feature. Again for purposes of illustration only, the discussion that 
follows will be couched in terms of traces contained within a stacked 3-D survey, although 
any assembled group of spatially related stacked or unstacked seismic traces could 
conceivably be used. 

Theoretica! Badcground 

In Figure 3 is illustrated a time series 310 together with various Fourier transform 
amplitude spectra computed therefrom. Spectrum 320 is computed over a relatively long 
window of seismic data (360), so its overall shape will tend to resemble die underlying 
seismic source wavelet. However, shorter windows 370 will tend yield spectra that are more 
affected by — and representative of — the underlying geology. (For a fiiUer discussion of 
this phenomena, see the invention taught by Gridley and Partyka in U.S. Patent No. 
5,870,691, the disclosure of which is incorporated haein by reference.) In either case, it 
should be clear to those of ordinary skill in the art that the calculated spectrum can take a 
variety of different shapes and those shapes might be a function of reflection time, position 
the surfece of the eardi, length of the analysis window, or many otiier variables. Although 
eitho- sort of analysis window (long or short) might used with the instant invention, in the 
preferred embodiment short windows will typically be used. 

Technical Discussion 

According to a first aspect of the instant invention, there is provided a method of 
seismic attribute generation for use in subsurface geophysical exploration which is based 
an analysis of the general shape of the frequency spectmm calculated from at least a potion of 
a seismic trace. 
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Figure 4 contains a broad overview of the principle steps in the preferred embodiment. 
As a first step, a seismic data set 410 containing seismic traces is accessed. This data set 
might be a stacked 2D line or 3D volume, a collection of unstacked traces, or any other 
arrangement of stacked or unstacked traces that are relatable to a position on the surface of 

5 the earth. For purposes of illustration hereinafter, it will be assumed that a stacked 3-D 
seismic volume is used in the analysis. 

The interpreter typically specifies a zone of interest within a particular 3-D volume. 
The zone of interest might be, by way of example, the undulating region bounded by two 
picked reflectors as is pictured in Figure 4. In that case, the reflector is preferentially 

10 flattened or datumized (i.e., made flat by shifting individual traces up or down in time, step 
413) before analysis, and possibly also palinspastically reconstructed. More conventionally, a 
specific bounded flat time interval (for example, fi-om 2200 ms to 2400 ms) might be 
specified, thereby defining a "cube" or, more accurately, a "box" of seismic data within the 3- 
D volume; a sub- volume. Note that the zone of interest "sub-volume" might be so expansive 

15 as to include the entire seismic volume. Additionally, the lateral extent of the zone of interest 
may be delimited by specifying bounding "in-lme" and "cross-line" trace limits. Other 
methods of specifying the zone of interest are certainly possible and have been contemplated 
by the inventor. The step of flattening the zone of interest 413 is not strictly required and 
would, of course, be unnecessary if the zone of interest were a "flat" slab of data. However, 

20 flattening an undulating zone of interest is a computational convenience for purposes of the 
steps that follow. 

As a next step, a seismic trace is selected 415 and the zone of interest within this trace 
is identified. Then the location of the analysis window is determined 420 within the zone of 
interest and the particular samples that will be used in the next step are identified. The 
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samples within the analysis window are then transformed via a Fourier transfomi and the 
amplitude spectrum 425 is computed (step 423) from the resulting transform coefiBcients. 
Note that a preferred variant of this step involves the computation of a composite amplitude 
spectrum which might be, by way of example only, the average of several short-window 
Fourier transform spectra that have been calculated within the same analj^is window or zone 
of interest. 

Additionally, it is preferable that a weightmg function should be applied to the 
samples in the analysis window before they are transformed. In flie preferred embodiment, 
the weight function wUl be a "Gaussian" weight One popular definition of the Gaussian 
weight fimction is: 

x2 . 



Wl 



(m) = or3e-('"-^> ^"2.m=0,L-l 
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where "L" is the length of the data window. The previous equation was given only as ; 
example and many other formulations are possible. In feet, any weight function that helps to 
mitigate the effects of the data window on the transform coeflRcients would be suitable. The 
previous steps are all well known to those skilled in the art and need not be discussed in any 
great detail herein. 

Given the calculated spectrum 425, as a next step a function characterized by constant 
coefficients is fit (steps 430 and 432) to the spectral values 430. In terms of equations, let the 
digital samples within the selected seismic trace be represented by the variable x(i), where "i" 
could potentially vary between 1 and the number of samples in the seismic trace, "N" 
hereinafler. However, for purposes of specificity, it will be assumed hereinafter that the 
samples are drawn fmm an analysis window that is shorter than the length of the seismic 
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trace, with the analysis window beginning at Mt and being "L" samples in length. The 
Fourier transform coefficients calculated from these values will be represented hereinafter as 
X(i), i = 0, \J1, and the amplitude spectrum A(i) is usually calculated as the complex 
magnitude of the corresponding Fourier transform coefficient: 

5 A(i) = lX(i)|,i = 0,L/2. 

As a next step, a functional form is selected for fitting to the amplitude spectrum. The 
preferred functional form is a low-order polynomial, but many other functions could be used 
instead. In terms of equations, the selected function, F(«), hereinafter, will have one or more 
constant coefficients ocq, a i, . . a m> which are to be estimated from the spectral data. The 

10 constants are to be selected such that 

A(f) = F(f;ao,ai, ...,aM) 
fits the spectral values, "A", as closely as possible, where inclusion of the parameter "f 
indicates that F(») is a function of spectral frequency. For example, the polynomial equation: 

A(f) = F(f;cto,ai, aM)= ao + aif++a2f^ 
15 describes a second order polynomial fit to the Fourier transform coefficients. In terms of 
matrices, the goal is to find a set of alpha coefficients that, as nearly as possible, fit the 
observed spectral amplitudes 



"A(l)" 




1 f(l) 


A(2) 




1 f(2) 


A(L). 




.1 f(L) 



|La2J 



where the quantity f(k) is preferably the "Icth" Fourier fiequency, 

k 



20 



f(k) = 



L( At/1000)' 



k = 0,f 
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and At is the sample rate in milliseconds. Alternatively, f(k) could be taken to be equal to the 
index "k" or any multiple thereof, since the Fourier frequencies are equally spaced for 
transforms calculated by a conventional discrete Fourier transform. Further, it might be 
desirable in some cases to transform the frequency axis to improve the quality the resulting 
fit. One way of doing that would be to set 

f(k) = kP 

where the exponent "p'' is chosen to in such as way as to expand or compress the frequency 
axis to improve the fit between the observed amplitudes and the chosen fimction. Nonlinear 
data transformations such as this — and others such as logarithms, powers, roots, etc. — are 
well known to those skilled in the art and need not be discussed further here. For fiirther 
examples of these sorts of transformations see, for example, Daniel and Wood, Fitting 
Equations to Data, Wiley & Sons, 1971, or Mosteller and Tukey, Data Analysis and 
Regression, Addison- Wesley, 1977, the disclosures of which are incorporated herein by 
reference. Finally, it should be noted that the transform frequencies f(k) need not be equally 
spaced, in which case the previous approach would need to be modified slightly to use the 
actual frequency values (rather than the index "k"). 

In terms of equations, the previous matrix expression may be written more compactly 

as: 

A«Fa, 

where tiie symbol has been used to indicate that the unknown constants (alpha) are to be 
chosen such that the left and right sides of the equation are as nearly equal as possible; where 
ttie matrix of Fourier frequencies (or transformations, thereof) has been represented as "F'; 
and where the vector of unknown coefficients has been designated as "a". It is well known to 
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those skilled in the art that under standard least squares theory, the choice of the coefficient 
vector that minimizes the difference between A and Fa is 

a = (A'^A + 8l)'^AF 

where the superscript indicates that a matrix inverse is to be taken. In the previous 
5 equation, a is a vector of coefiBcient estimates, the quantity I represents the M by M identity 
matrix (3 by 3 in the previous example), and e is a small positive number which has been 
introduced — as is commonly done — for purposes of stabilizing the matrix inversion. 
Finally, those skilled in the art will recognize that the least squares minimization of the 
previous matrix equation is just one of many nonns that might be used to constrain the instant 
10 problem and thereby yield a solution in terms of the unknown alphas, some alternatives 

norms being, by way of example only, the Li or least absolute deviation norm, the Lp or least 

"p" power norm, and many other hybrid norms such as those suggested in the statistical 
literature on robust estimators. See, for example, Peter J. Huber, Robust Statistics^ Wiley, 
1981, the disclosure of which is incorporated herein by reference. Further, non-linear 
1 5 equations such as : 

may also be used, although the previous (linear) method of solving this sort of (nonlinear) 
equation would need to be modified according to methods well known to those skilled in the 
art, 

20 Further, it should be noted that instead of using the spectral amplitude of each 

transformed analysis window, any other quantity that is computed ifrom the Fourier transform 
coefficients could have been used instead. For example, a function could be fit to the phase 
spectrum (preferably after first "unwrapping the phase" as that term is known to those skilled 
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in the art). Alternatively, the original complex Fourier transform coefficients could have been 
used in the fit, although the previous matrix might need to be separated into separate "real" 
and "imaginary" equations in order to obtain a solution. Still fiirther, smoothed spectra 
(including the integrated spectra or a cumulative sum of the spectral values) could also be 
used as could mathematical functions of the individual spectral or Fouri^ transform 
coefficient values (e.g., square roots, logarithms, inverse, roots, numerical derivatives, etc.). 
In the preferred embodiment, the cumulative amplitude spectrum is used in the fitting process 
and the preferred fitting function is a first order polynomial (i.e., a linear equation). 

The alpha coefficient estimates that are produced by this process provide a wealth of 
seismic attributes which may be mapped and analyzed by the explorationist: at least as many 
attributes are produced as there are constants m the fitted equation. La the preferred 
embodiment, all of the coefficients that are calculated within each analysis window are 
written to output as separate planes (in the case of 3D data) of output. Thus, if a linear 
equation is fit to the spectra, at least two planes of coefficients would preferably be produced. 
Additionally, it is preferred that the residual to the fit (e.g., the parameter "SSE" as discussed 
below) is also written to output, thus yielding a third attribute "slice" that is preferably 
produced from a fit of a linear functional form. 

Some examples are given hereinafter of how these coefficients may be used in 
practice, but 

the suggestions detailed below represent only a few of the many uses to which these 
coefficients may be put. In the event that a second order polynomial is fit to the amplitude 
spectrum, the sign of the coefficient associated with the quadratic term (a2 in the previous 
example) will tend to reflect whether the spectrum is "bowl" shaped (e.g., with a pronounced 
notch in the spectrum, e.g.. Figure 3, item 330) or "hump" shaped (as it is illustrated in Figure 

19 



wo 01/96905 



PCT/USOl/40506 



4, item 425), with a positive estimated coefficient value tending to indicate a spectrum having 
a bowl shape. Similarly, if a first order (linear) equation is fit to the unwrapped phase, the 
linear regression component of the fit will tend to reflected any time shifts that are present in 
the data, it being well known that a time-domain shift is expressed in the phase spectrum as a 

5 linear trend. Clearly, it would be impractical to list here all of the possible interpretations of 
the many different calculated coefficients (and mathematical combinations of those 
coeflBcients) that might be used: the interpretation of each coefficient will necessarily be 
based on the particular function fitted. Those of ordinary skill in the art will be able to device 
any number of fitting functions and corresponding coefficient interpretations, 

10 Another quantity that is may be calculated in connection with the fitting process is the 

degree to which the data are adequately *'fit" by the selected equation within each analysis 

2 

window. The statistical Coefficient of Determination, "r" or "r " as that quantity is known to 

those skilled in the art, is one conventional measure of the "goodness of fit" of selected 

2 

equation to the data (e.g., spectral values). A large value of r (i.e., a value near 1.0) indicates 

15 that the data are well fit withm that window level by the selected equation. However, a value 

of r near zero signifies that the Fourier transform coefficients analyzed within that particular 

window are not well fit by the selected equation. Statistical coefficients related to the quality 
of the fit will be considered hereinafter to be another fitted constant parameter. 

According to anoAer aspect of the instant invention, there is provided a method 
20 substantially as described above, but wherein display of the calculated coefficient is 

contingent on the value of the correlation or on some other goodness of fit parameter. In the 
preferred embodiment, the seismic data will be pre-processed as described previously and a 
zone of interest and analysis window length will be selected by the seismic processor. Then a 
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trace is read, the data samples in the analysis window are extracted, and the samples are 
transformed via a Fourier transform as before. The spectrum is then formed and is fit using 
the functional form of choice and, additionally, the coefficient of determination is calculated 
for the fit. Howver, rather than directly storing or displaying the calculated coefficient as a 
seismic attribute, the Coefficient of Detennination is firet compared with some predetemiined 
value. Then, dependmg on the results of the comparison, either the calculated value or some 
other value will be displayed. For example, consider the following rule that displays the 
original coefficient estimate if tiie Coefficient of Detennination is "high" (i.e., the fit is 
relatively good) and the value "zero" if the coefficient is "low": 



[0,else 



y(Mt)= 

where y(Mt) is a storage area (typically an empty seismic trace) for receiving the seismic 
attribute calculated from an analysis window that begins at sample Mt; where ? is the 
Coefficient of Determination for the fit for of a spectrum calculated fo)m an analysis window 
that begins at time M, , where az is the calculated estimate of the parameter 0.2, and where "d" 
is an arbitrary positive value selected in acco«lance with the discussion that follows. 

Those skilled in the art know that in regression analysis, the Coefficient of 
Determination is a statistical measure of causal association between two variables, where the 
higher the numerical value (on a scale of 0.0 to +1.0) the greater the degree of association. 
More formally, this statistic is a measure of the proportion of the total variation in the 
predicted variable that is associated with variations in tiie predictore in the regression 
equation. Thus, a Coefficient of Determination o^ say, 0.85 means that 85% of the variance 
in the predicted variable is attributable to tiie regression model, whereas 15% is unexplained. 
In brief, values of tiie Coefficient of Determination near 1.0 indicate that tiie data are well fit 
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by the regression model, whereas values near zero indicate that a poor fit with the model. 
Although there are many ways to express the computational formula for the Coefficient of 
Determination, a preferred way to do so is as follows: 



5 where X is the average of the spectral values X(*) within an analysis window of length "L"; 

where ao, ai, and a2 are estimated alpha coefficient values (for the quadratic example 
discussed previously); where "Mean SS" is the sum of squares about the mean (the first term 
in the numerator); and, "SSE" is the sum of squares due to lack of fit by the regression (the 
second term in the numerator). Those skilled in the art will recognize that there are many 

10 other ways to formulate the Coefficient of Determination. However, the essential feature of 
this statistical measure is that it returns some indication of the goodness of fit of the chosen 
fiinctional form to the spectral data within that analysis window. 

In terms of the instant embodiment, relatively larges values of the Coefficient of 
Determination indicate that the seismic data at those time levels are well fit by the chosen 

15 fiinctional form. Thus, in this instant embodiment, the calculated coefficient is presented to 
the explorationist only for those analysis windows tthat have an associated r^ that is high, i.e., 
"near" 1.0. One preferred choice for the parameter "d" is about 0.50, although the value that 
works best for a given functional form and data set may need to be determined by trial and 
error. In the alternative, of course, the calculated coefficient might be presented to the 




m=0 



MeanSS-SSE 
Total SS 
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e3q)loratiomsEt only for those analysis windows with values of near zero, however the 
interpretation of tfie resulting display would then be different. 

For those coefiBcients that correspond to low values of r^ a constant replacement 
value is preferably used to indicate that feet, h the preferred embodiment, and as is 
illustrated in the previous equation, the value zero has been selected as the substitute value. 
However, "zero" was selected for purposes of convenience only and many other choices 
might be used instead. A main purpose of using a single constant replacement value is to 
make inferior functional fits readily recognizable by the explorationist who will be viewing 
and analyzing the output from this method. 

A display of seismic attributes formed in this manner provides a means for the 
explorationist to rapidly identify within a section or volume those locations in time and space 
where fte data are accurately fit by the selected fiinctional form. ITiis display would also 
allow the explorationist to quickly locate regions of poor fit throughout an entire seismic 
section or volume. 

OTHER BASIS FUNCTIONS AND WEIGHT FUNCTIONS 
Those of ordinary skill in the art will recognize that, although the present invention is 
discussed herein in terms of the discrete Fourier transform, in reality the Fourier transform is 
just one ofmiy number ofdiscrete time data transformations that could used in exactly*^ * . 
same fashion. The general steps of (1) selecting an analysis window; (2) calculating a data 
transformation from the samples within that window; (3) fitting a fimction characterized by 
constant coefiBcients to the resulting transfonn coeflFicients (or mathematical fiinctions of 
those coefficients); and, (4) selecting one of the coefficient estimates to use as a seismic 
attribute, could be accomplished with a wide variety of discrete data transformations other 
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than the Fourier. If the transformation is other than a Fourier, the only real changes are that 
the resulting transform coefficients / spectrum will be interpreted in terms of the underlying 
basis functions (rather than in terms of frequency) and changes in the shape of transform 
coefficients may be also interpreted differently. 

5 Those skilled in the art will understand that a discrete Fourier transform is just one of 

many discrete linear unitary transformations that satisfy the following properties: (1) they are 
linear operators that are (2) exactly invertible, and (3) their basis functions form an 
orthonormal set. For the discrete Fourier transform, the basis functions corresponding to a 
forward transform of length L are conventionally chosen to be the set of complex 

10 exponentials and each transform coefficient corresponds to the strength of a particular a 

frequency in the data. By way of a specific example of how an alternative transform could be 
used, those skilled in the art understand that a discrete Walsh transform could be used in 
place of the Fourier transform and the Walsh coefficients similarly fitted, displayed, and 
analyzed. In the manner disclosed above, a Walsh transform may be computed within an 

15 overlapping series of sliding windows and the coefficients resulting therefrom subjected to a 
functional fit as described previously. Rather than the calculated transform coefficients 
representing frequency, of course, these coefficients instead represent a similar quantity called 
"sequency" by those skilled in the art. Other examples of transforms that might prove to 
useful in particular settings include, without limitation. Hartley transforms, Fourier Sine and 

20 Fourier Cosine transforms, and wavelet transforms. 

Finally, the weight / tapering function discussed previously could take any number of 
forms. Some of the more popular data windows are the Hamming, Hanning, Parzen, Bartlett, 
and Blackman widows. Each of these functions has certain advantages and disadvantages. 
However, the Gaussian taper is in many ways optimal for this application. 
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Conclusions 

From the foregoing description, it will be observed that numerous variations, 
alternatives and modifications will be apparent to those skilled in the art. Accordingly, this 
description is to be construed as illustrative only and is for the purpose of teaching those 
skilled in the art the manner of cariying out the invention. By way of example, other 
algorithms may be used to calculate the right pseudo-inverse matrix. Moreover, equivalent 
computations may be substituted for those illustrated and described. Also certain features of 
the invention may be used independently of other features of the invention. 

In the previous discussion, the language has been expressed in terms of operations 
performed on conventional seismic data. But, it is understood by those skilled in the art that 
the invention herein described could be applied advantageously in oflier subject matter areas, 
and used to locate oflier subsurface minerals besides hydrocarbons. By way of some 
examples of trace data other than conventional seismic data that would be suitable for use 
with tfie invention disclosed herein, the same approach described herein could be used to 
process and/or analyze multi-component seismic data, shear wave data, converted mode data, 
VSP data, magneto-telluric data, cross well survey data, full waveform sonic logs, or model- 
based digital simulations of any of the foregoing. Additionally, the mefliods claimed 
hereinafter can be applied to mathematically transformed versions of these same data traces 
including, for example: frequency domain Fourier transformed data; transformations by 
discrete orthonormal transforms; wavelet-transformed data; instantaneous phase traces, 
instantaneous frequency traces, the analytic and quadrature traces; etc. Additionally, tfie 
methods claimed herein after can be applied to mathematically transformed versions of these 
same data traces including, for example: frequency domain Fourier transformed data; 
transformations by discrete orthonormal transforms; instantaneous phase traces, instantaneous 
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frequency traces, quadrature traces, and analytic traces, etc. In short, the process disclosed 
herein can potentially be applied to any geophysical time series, but it is preferably applied to 
a collection of spatially related time series containing structural or stratigraphic information. 
Thus, in the text that follows those skilled in the art will understand that "seismic trace" is 

5 used herein in a generic sense to apply to geophysical time series in general. 

While the inventive device has been described and illustrated herein by reference to 
certain preferred embodiments in relation to the drawings attached hereto, various changes 
and fiirther modifications, apart from those shown or suggested herein, may be made therein 
by those skilled in the art, without departing from the spirit of the inventive concept, the 

10 scope of which is to be determined by the following claims. 
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THE CLAIMS: 

WHAT IS CLAIMED IS: 

1 . A method of seismic exploration, 

wherein there is provided a seismic survey consisting of seismic traces 
containing digital samples collected over a predetermined volume of tfie earth 
containing subsurface features conducive to the generation, migration, 
accumulation, or presence of hydrocarbons, 
comprising the steps of: 

(a) selecting a particular seismic trace, said seismic trace havmg at least two 
digital samples contained therein; 

(b) choosing an analysis window containing at least two of said at least two digital 
samples; 

(c) transforming said at least two digital samples using a discrete orthonormal 
transformation, said discrete orthonormal transformation producing transform 
coefficients from said digital samples so transformed; 

(d) fitting a function characterized by at least one constant coefficient to said 
transform coefficients, thereby producing at least one coefficient estimate; and, 

(e) forming a seismic attribute from any from any coefficient estimates so 
calculated, said seismic attribute for using in identifying said subsurface 
features conducive the generation, migration, accumulation, or presence of 
hydrocarbons. 



2. A metfiod of seismic exploration according to Claim 1, wherein step (d) includes the 
steps of: 
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(dl) calculating an amplitude spectrum from said transform coefficients, 
and, 

(d2) fitting a function characterized by at least one constant coefficient to 
said amplitude spectrum, thereby producing at least one coeflBcient 
estimate. 

3. A method of seismic exploration according to Claim 1, wherein step (d) includes the 
steps of: 

(dl) calculating a phase spectrum from said transform coefficients, and, 
(d2) fitting a function characterized by at least one constant coefficient to 

said phase spectrum, thereby producing at least one coefficient 

estimate. 

4. A method of seismic exploration according to Claim 1, further including the step of: 
(f) storing said seismic attribute. 

5. A method of seismic exploration according to Claim 4, wherein steps (a) through {f) 
are repeated a predetermined number of times. 

6. A mediod of seismic exploration, wherein are provided stored seismic attributes 
prepared according to the method of Claim 5, comprising the steps of: 

(a) accessing said stored seismic attributes; and, 

(b) displaying said stored seismic attributes. 
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7. A method of seismic exploration, 

wherein there is provided a seismic survey consisting of seismic traces 
containing digital samples collected over a predetermined volume of the earth 
containing subsur&ce features conducive to the generation, migration, 
accumulation, or presence of hydrocarbons, 
comprising the steps of: 

(a) selecting a particular seismic trace, said seismic trace having at least two 
digital samples contained therein; 

(b) choosing an analysis window containing at least two of said at least two digital 
samples; 

(c) transforming said at least two digital samples using a discrete orthonormal 
transformation, said discrete orthonormal transformation producing transform 
coefiBcients from said digital samples so transformed; 

(d) fitting a function characterized by at least one constant coefficient to said 
transform coefiBcients; 

(e) determining at least one goodness of fit value corresponding to said functional 
fit of step (d); and, 

forming a seismic attribute fi*om any from any goodness of fit values so 
determined, said seismic ath-ibute for using in identifying said subsurface 
features conducive the generation, migration, accumulation, or presence of 
hydrocarbons. 
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